!--------------------------------------------------------------------------------------
subroutine vmd_init_modulation(xvec, kXvec, xMod, uPerpMod)
!--------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward
  use dimensions, only : jx, nkx
  use parameters, only : Lx, SpeckleWidth 

  implicit none

  real,    intent(in)  :: kXvec(jx)
  real,    intent(in)  :: xvec(jx)
  complex, intent(out) :: xMod(jx), uPerpMod(jx)


!-- local variables ---

  integer              :: ix, ikx
  

!-- create modulation vectors --------------------------------------------------------


      xMod(:) = exp( - 20 * xvec(:)**2 / SpeckleWidth**2 )   !  complex xMod(jx)
 
      call fft_forward(jx,xMod,xMod)
      xMod=xMod/float(jx)
      call deAlias(jx,xMod)
 
      open(3,file='phi0_vs_kx.txt')
      write(3,*) '% kx <|phi0(kx)|^2>'
      do ikx=1, nkx/2  - nkx/6 -1 ; 
         write(3,501) kxVec(ikx), abs( xMod(ikx) )**2
      end do
      close(2)


      call fft_backward(jx,xMod,xMod)  

      open(26,file='Potential_xMod.txt'); 
      write(26,*) '%x    xMod'
      do ix =1, jx
         write(26,501) xvec(ix), real(xMod(ix))
      enddo
      close(26)
        
 
      uPerpMod(:) = exp( -14 * xvec(:)**6 / (Lx/2)**6 )   !  complex uPerpMod(jx)
 
      call fft_forward(jx,uPerpMod,uPerpMod)
      uPerpMod=uPerpMod/float(jx)
      call deAlias(jx,uPerpMod)        
      call fft_backward(jx,uPerpMod,uPerpMod)
        

      open(5, file='output', position='append')
      write(5,*) 'uPerpMod(1 ) =', uPerpMod(1 )
      close(5)

      open(26,file='uPerpMod.txt')
      write(26,*)'%x    uPerpMod'
      do ix =1, jx 
         write(26,501) xvec(ix), real(uPerpMod(ix))
      enddo
      close(26)

501   format(  6400(e12.6,2x)  )


end subroutine vmd_init_modulation

!--------------------------------------------------------------------------------------
subroutine vmd_init_density(zvec, ExFilter, xmod, uPerpMod, E0C,  &  ! inputs
                            E0z, phi02DFT)                           ! outputs
!--------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward
  use dimensions, only : jx, jz
  use parameters, only : phi0, kz, iE0Modulation, iunit


  implicit none

  real,    intent(in)  :: zvec(jz)
  real,    intent(in)  :: ExFilter(jx)
  complex, intent(in)  :: xMod(jx)
  complex, intent(in)  :: uPerpMod(jx)
  complex, intent(in)  :: E0C(jz)

  real,    intent(out) :: E0z(jx,jz)
  complex, intent(out) :: phi02DFT(jx,jz)

 
!-- local variables ---

  integer              :: ix, iz       ! counters in r-space
  integer              :: ikx           ! counters in k-space

  complex              :: Phi1Dx(jx)
  complex              :: density1D(jz)
  complex              :: phi02D(jx,jz)

!-- modulate source -------------------------------------------------------------------
  

      !  modulate INITIAL source in x direction
      !  E0C(iz) =  -iunit * kz * phi0 *exp ( iunit * kz * (iz-1) * dz  )     
      !  so when x=SpeckleWidth/2, E0z is down by a factor of 0.00674
      !  E0C(iz) = -iunit * kz * phi0 *exp ( iunit * kz * z  )

      if (iE0Modulation==1) then
         do ix=1,jx
           E0z(ix,:)    = xMod(ix) * E0C(:) 
           phi02D(ix,:) = xMod(ix) * phi0 *exp ( iunit * kz * zvec(:)  ) 
        enddo
      else
         do ix=1,jx
           E0z(ix,:)    = E0C(:)
           phi02D(ix,:) = phi0 *exp ( iunit * kz * zvec(:)  )
         enddo
      endif  


      do iz=1,jz
         Phi1Dx(:) = phi02D(:,iz)       !    density2D(all x,kz)   complex Phi1Dx(jx)
         call fft_forward(jx,Phi1Dx,Phi1Dx)
         Phi1Dx=Phi1Dx/float(jx)
         call deAlias(jx,Phi1Dx)
         phi02DFT(:,iz) = Phi1Dx(:) * ExFilter(:) !   ph02DFT( all kx, all z)
      enddo

      do ikx=1,jx
         density1D(:) = phi02DFT(ikx,:)
         call fft_forward(jz,density1D,density1D)
         density1D=density1D/float(jz)
         call dealias(jz,density1D)
         phi02DFT(ikx,:) = density1D(:)   !   ph02DFT( all kx, all kz)
      enddo

!-------------------------------------------------------------------------------------

end subroutine vmd_init_density
